##############
# Results.R

# This document creates all plots and tables for the 
# Results section of the Austrian election Study

# Library
library(ggplot2)
require(MASS)
library(xtable)
library(reshape)
library(reshape2)
library(rjags)

# Load Data
source("Model_Data_oest.R")
load("oest_long.Rdata")
load("models_oest_rat.Rdata")


##############################
# Figure 2

cq <- function(res, var = "gamma" ){
  
 res2 <-   as.matrix(res)
 sel <- grep(var,colnames(res2))
 ret <-  t(apply(res2[,sel],2,quantile,c(0.05,0.5,0.95)))
  
 rownames(ret) <- rep("",nrow(ret))
 return(ret)
 
  }

# Get First-difference 

  df <- data.frame(cbind(as.matrix(ovp.gruene)[,"diff"]
                         ,as.matrix(spo.gruene)[,"diff"]
                         ,as.matrix(ovp.fpo)[,"diff"]
                         ,as.matrix(spo.fpo)[,"diff"]
  ))
  
  
  df <- as.data.frame(t(apply(df,2,quantile,c(0.01,0.5,0.99))))
  colnames(df) <- c("low","mid","high")
  
  df$coal <- as.factor(1:4)
  levels(df$coal) <- c("ÖVP-Greens","SPÖ-Greens","ÖVP-FPÖ","SPÖ-FPÖ")
  
  
  df$type <- 2
  df$choice <- 3


  # Get Weights 
  df2 <- data.frame(rbind(cq(ovp.gruene)
                         ,cq(spo.gruene)
                         ,cq(ovp.fpo)
                         ,cq(spo.fpo)))
  
  
  names(df2) <- c("low","mid","high")
  df2$coal <- as.factor(sort(rep(1:4,2)))
  levels(df2$coal) <- c("ÖVP-Greens","SPÖ-Greens","ÖVP-FPÖ","SPÖ-FPÖ")

  
  df2$type <- 1
  df2$choice <- rep(1:2,2)

  # Combine Both
  
  df <- rbind(df,df2)

  df$type <-  as.factor(df$type)
  levels(df$type) <- c("Mixing parameters","First differences") 
  
  df$choice <-  as.factor(df$choice)
  levels(df$choice) <- c("Standard","Vignette","First difference") 
  
  
  # Scales

  pdf("Oest-plot-joint.pdf",height=4)
  ggplot(df, aes(y=mid,ymin=low,ymax=high,x=1,col=choice)) + 
    geom_pointrange(size=1,position=position_dodge(0.5))  + coord_flip() + 
    facet_grid(coal~ type) +
    scale_color_grey() +  
    scale_x_continuous(name="", labels="" , breaks=1,limits=c(0,2)) + theme_bw() + 
    theme(legend.direction = "horizontal", legend.position = "top", legend.title=element_blank(), 
          strip.text.y = element_text(size=10,angle=360)) +
    scale_y_continuous(name="", breaks=c(0,0.25,0.5,0.75,1), 
                       labels=c("0","0.25","0.5","0.75","1"),limits=c(0,1)) 
  dev.off()

#############################
# Figure 3: Simulation Plot

controls <- c("PID","gender","educ","relig","union","income","age" )

sgd  <- Model.data(party = "spo" ,coalition = "gruene"
                  ,vig = "wahl2"
                  ,rat.coal="rat.spo.gruene"
                  ,covariates=controls
                  )


# Function simulates qunataties of interest for a given set of covariates

pf <- function(x,S1=S,a1=1,a2=0,X=cova,PID=P,vote=2,lr=LR,qunat=T,return_p1 = FALSE) {
  
  alpha  <- S1[,grep("gamma",colnames(S1))]
  delta <- S1[,grep("delta",colnames(S1))]
  beta <- cbind(S1[,"beta[1,1]"], S1[,"beta[1,2]"],S1[,"beta[2,1]"],S1[,"beta[2,2]"])
  psi <- S[,grep("psi",colnames(S))]
  
  sel.x11 <- grep("psi\\[1,1",colnames(psi))
  sel.x21 <- grep("psi\\[2,1",colnames(psi))
  sel.x12 <- grep("psi\\[1,2",colnames(psi))
  sel.x22 <- grep("psi\\[2,2",colnames(psi))
  
  U11 <- beta[,1] + delta[,1]* (alpha[,1]*x[1] + (1-alpha[,1]) *x[3])  +delta[,6]*lr[1] + delta[,3] * PID[1]  + psi[,sel.x11] %*% X
  U21 <- beta[,3] + delta[,2]* (alpha[,2]*x[1] + (1-alpha[,2]) *x[3]) +delta[,7]*lr[1] +  delta[,4] * PID[1]  + delta[,5] * a1  +psi[,sel.x21] %*% X
 
  U12 <- beta[,2] + delta[,1]* (alpha[,1]*x[2] + (1-alpha[,1]) *x[3])  +delta[,6]*lr[2]+ delta[,3] * PID[2] + psi[,sel.x12] %*% X
  U22 <- beta[,4] + delta[,2]* (alpha[,2]*x[2] + (1-alpha[,2]) *x[3])  +delta[,7]*lr[2]+ delta[,4] * PID[2] + delta[,5] * a2 + psi[,sel.x22] %*% X

  U13 <- 0
  U23 <- 0

  # P1
  sumU1 <- exp(U11)+exp(U12)+exp(U13)
  p12 <- exp(U12)/sumU1
  p11 <- exp(U11)/sumU1

  # P2
  sumU2 <- exp(U21)+exp(U22)+exp(U23)
  p22 <- exp(U22)/sumU2
  p21 <- exp(U21)/sumU2


  if(vote==2){
    # qunatiles 
    if(qunat){
      return(c(mean(p22),quantile(p22,c(0.025,0.975)),mean(p12),quantile(p12,c(0.025,0.975))))
    } else{
      if(return_p1) return(p12) else return(p22)
    }
  } 
  if(vote==1) {
    # qunatiles
    if(qunat){
      return(c(mean(p21),quantile(p21,c(0.025,0.975)),mean(p11),quantile(p11,c(0.025,0.975))))
    } else {
      if(return_p1) return(p11) else return(p21)
    }
  }
  # return(p22-p12) to caclulate drop in predicted probabilities
  }

 
 
# AVERAGE VALUE APPROACH: Set Covariates
cova <- apply(as.matrix(sgd[controls[-1]]),2,mean)
Rat <- cbind(7,11,seq(1,11,1))
P <- c(0,0)
LR <- apply(as.matrix(sgd[c("lr1","lr2")]),2,mean)
S <- as.matrix(spo.gruene[[1]])

# Get Draws for different scenarios
d <- data.frame(t(apply(Rat,1,pf,a1=0,a2=1,vote=2,qunat=T)))


# Put into Data-frame
colnames(d) <- c("mean.1","low.1","high.1","mean.2","low.2","high.2")
d$x <- seq(1,11,1)
d <- reshape(d,direction="long",varying=1:6)
d$Vote <- as.factor(d$time)
levels(d$Vote) <- c("Vignette","Standard")


pdf("Sim-plot.pdf",height=5)

ggplot(d,aes(x=x,y=mean,ymin=low,ymax=high,group=time,fill=Vote)) + ylab("Probability of voting for the Green party") + 
  geom_line(aes(linetype=Vote))  +  geom_ribbon(alpha=0.5,fill="grey")  + scale_colour_grey() +  theme_bw() + labs(colour = "")   + 
  xlab("Ratings for SPÖ-Greens coalition") + scale_x_continuous(breaks=c(1,11),labels=c("low","high"))   + 
  theme(legend.justification=c(1,0), legend.position=c(1,0))   

dev.off()



#############################
# Figure 5: Simulation Plot for all models


################
# OGD

ogd  <- Model.data(party = "ovp" ,coalition = "gruene"
                  ,vig = "wahl2"
                  ,rat.coal="rat.ovp.gruene"
                  ,covariates=controls
                  )

cova <- apply(as.matrix(ogd[controls[-1]]),2,mean)
LR <- apply(as.matrix(ogd[c("lr1","lr2")]),2,mean)
S <- as.matrix(ovp.gruene[[1]])

# OVP Voter
#############

Rat <- cbind(11,7,seq(1,11,1))
P <- c(1,0)


d <- data.frame(t(apply(Rat,1,pf,a1=1,a2=0,vote=1)))
colnames(d) <- c("mean.1","low.1","high.1","mean.2","low.2","high.2")
d$x <- seq(1,11,1)

d <- reshape(d,direction="long",varying=1:6)
d$Vote <- as.factor(d$time)
levels(d$Vote) <- c("Vignette","Standard")
d$Rat <- seq(1,11,1)
d$Vote <- "OVP"
d$Vig <- "OVP-Green"

# Green Voter
##################

Rat <- cbind(7,11,seq(1,11,1))
P <- c(0,1)
d2 <- data.frame(t(apply(Rat,1,pf,a1=0,a2=1)))
colnames(d2) <- c("mean.1","low.1","high.1","mean.2","low.2","high.2")
d2$x <- seq(1,11,1)

d2 <- reshape(d2,direction="long",varying=1:6)
d2$Vote <- as.factor(d2$time)
levels(d2$Vote) <- c("Vignette","Standard")



d2$Rat <- seq(1,11,1)
d2$Vote <- "Green"
d2$Vig <- "OVP-Green"

d <- rbind(d,d2)


################
# SGD


sgd  <- Model.data(party = "spo" ,coalition = "gruene"
                  ,vig = "wahl2"
                  ,rat.coal="rat.spo.gruene"
                  ,covariates=controls
                  )

cova <- apply(as.matrix(sgd[controls[-1]]),2,mean)
LR <- apply(as.matrix(sgd[c("lr1","lr2")]),2,mean)
S <- as.matrix(spo.gruene[[1]])

# SPO Voter
#############

Rat <- cbind(11,7,seq(1,11,1))
P <- c(1,0)


d2 <- data.frame(t(apply(Rat,1,pf,a1=1,a2=0,vote=1)))
colnames(d2) <- c("mean.1","low.1","high.1","mean.2","low.2","high.2")
d2$x <- seq(1,11,1)

d2 <- reshape(d2,direction="long",varying=1:6)
d2$Vote <- as.factor(d2$time)
levels(d2$Vote) <- c("Vignette","Standard")

d2$Rat <- seq(1,11,1)
d2$Vote <- "SPO"
d2$Vig <- "SPO-Green"


d <- rbind(d,d2)


# Green Voter
##################

Rat <- cbind(7,11,seq(1,11,1))
P <- c(0,1)
d2 <- data.frame(t(apply(Rat,1,pf,a1=0,a2=1)))
colnames(d2) <- c("mean.1","low.1","high.1","mean.2","low.2","high.2")
d2$x <- seq(1,11,1)

d2 <- reshape(d2,direction="long",varying=1:6)
d2$Vote <- as.factor(d2$time)
levels(d2$Vote) <- c("Vignette","Standard")

d2$Rat <- seq(1,11,1)
d2$Vote <- "Green"
d2$Vig <- "SPO-Green"

d <- rbind(d,d2)



################
# OFD

ofd  <- Model.data(party = "ovp" ,coalition = "fpo"
                  ,vig = "wahl2"
                  ,rat.coal="rat.ovp.fpo"
                  ,covariates=controls
                  )

cova <- apply(as.matrix(ofd[controls[-1]]),2,mean)
LR <- apply(as.matrix(ofd[c("lr1","lr2")]),2,mean)
S <- as.matrix(ovp.fpo[[1]])

# OVP Voter
#############

Rat <- cbind(11,7,seq(1,11,1))
P <- c(1,0)
d2 <- data.frame(t(apply(Rat,1,pf,a1=1,a2=0,vote=1)))
colnames(d2) <- c("mean.1","low.1","high.1","mean.2","low.2","high.2")
d2$x <- seq(1,11,1)

d2 <- reshape(d2,direction="long",varying=1:6)
d2$Vote <- as.factor(d2$time)
levels(d2$Vote) <- c("Vignette","Standard")

d2$Rat <- seq(1,11,1)
d2$Vote <- "OVP"
d2$Vig <- "OVP-FPO"

d <- rbind(d,d2)

# FPO Voter
##################

Rat <- cbind(7,11,seq(1,11,1))
P <- c(0,1)
d2 <- data.frame(t(apply(Rat,1,pf,a1=0,a2=1)))
colnames(d2) <- c("mean.1","low.1","high.1","mean.2","low.2","high.2")
d2$x <- seq(1,11,1)

d2 <- reshape(d2,direction="long",varying=1:6)
d2$Vote <- as.factor(d2$time)
levels(d2$Vote) <- c("Vignette","Standard")

d2$Rat <- seq(1,11,1)
d2$Vote <- "FPO"
d2$Vig <- "OVP-FPO"

d  <- rbind(d,d2)

################
# SFD


 sfd  <- Model.data(party = "spo" ,coalition = "fpo"
                   ,vig = "wahl5"
                   ,rat.coal="rat.spo.fpo"
                   ,covariates=controls    
                   )

cova <- apply(as.matrix(sfd[controls[-1]]),2,mean)
LR <- apply(as.matrix(sfd[c("lr1","lr2")]),2,mean)
S <- as.matrix(spo.fpo[[1]])

# SPO Voter
#############

Rat <- cbind(11,7,seq(1,11,1))
P <- c(1,0)
d2 <- data.frame(t(apply(Rat,1,pf,a1=1,a2=0,vote=1)))
colnames(d2) <- c("mean.1","low.1","high.1","mean.2","low.2","high.2")
d2$x <- seq(1,11,1)

d2 <- reshape(d2,direction="long",varying=1:6)
d2$Vote <- as.factor(d2$time)
levels(d2$Vote) <- c("Vignette","Standard")

d2$Rat <- seq(1,11,1)
d2$Vote <- "SPO"
d2$Vig <- "SPO-FPO"


d <- rbind(d,d2)

# FPO Voter
##################

Rat <- cbind(7,11,seq(1,11,1))
P <- c(0,1)
d2 <- data.frame(t(apply(Rat,1,pf,a1=0,a2=1)))
colnames(d2) <- c("mean.1","low.1","high.1","mean.2","low.2","high.2")
d2$x <- seq(1,11,1)

d2 <- reshape(d2,direction="long",varying=1:6)
d2$Vote <- as.factor(d2$time)
levels(d2$Vote) <- c("Vignette","Standard")

d2$Rat <- seq(1,11,1)
d2$Vote <- "FPO"
d2$Vig <- "SPO-FPO"

d  <- rbind(d,d2)


#######################
# Graph

d$Vote <- as.factor(d$Vote)
d$Vig <- as.factor(d$Vig)
d$Choice <- as.factor(d$time)
levels(d$Choice) <- c("Vignette","Standard") 

pdf("Sim-plot3.pdf",height=5)
ggplot(d,aes(x=Rat,y=mean,ymin=low,ymax=high,group=Choice)) + ylab("Probability of voting for .... ") + 
  geom_line(aes(linetype=Choice))  +  geom_ribbon(alpha=0.2)  +  theme_bw() + 
  labs(colour = "")   + xlab("Coalition Ratings") + 
  scale_x_continuous(breaks=c(1,11),labels=c("low","high"))   + 
  theme(legend.justification=c(1,0), legend.position=c(0.2,0))   + 
  facet_grid(Vig  ~ Vote)
dev.off()

